% log-likehood function
% required program: T2_fun3, cpdf_wv_integrand3

function [out, lowu] = ll_fun5(params,data,I)
% params must be COLUMN vector with the parameters ordered as follows:
% (beta, lambda, muu, sigma2u) = params(1:4)
% alphav = params(5:ng+4)
% betav = params(I+5:2*ng+4)
% (alphaw1, alphaw2) = params(2*ng+5:2*ng+6)
% barw = params(2*ng+7)

% data = [dataP dataQ] with dataP and dataQ of size Tx1 and TxI, resp.

% muv0 must be ROW vector with initial values for muV

% ng is scalar - number of groups. currently 6. 
% I is scalar (number of firms/countries)
% EXAMPLES:
% params = [beta_d lambda_cost meanU sigmaU^2 alphaV betaV 1 1]
% ll_fun5(params',data,muV,I)

temp0c = params(2)+(I+1)*params(1);
tempD = 1/(params(2)/params(1) + I+1);

% T1 = [sum(params(1:2))/temp0c params(1)*ones(1,I); ones(I,1)/temp0c -eye(I)];
T1inv = [1 tempD*ones(1,I)/params(1); temp0c*tempD*ones(I,1) -eye(I)+tempD*ones(1,1)]';
%inv_absdetT1 = abs(1 + [0 params(1)*ones(1,I)*temp0c/sum(params(1:2))]...
%    *T1inv*[0; ones(I,1)/temp0c])*temp0c/sum(params(1:2));  %%CHECK THIS. WHY INV? DET(T1)=1?

T = size(data,1);
utildev_fun = @(t) (data - repmat(T2_fun3(params(1),params(2),t,I),T,1))*T1inv';
tempM1 = eye(I+1);  % tempM2 = tempM1(:,2:I+1);

% choice of w
barw = params(2*I+7);
int_correc = 0.025;    % to avoid Inf integrals
int_lim = barw/temp0c;
ntgrid = 25; htgrid = 2*int_lim/ntgrid; % grid for integration; check later
tgrid = (-int_lim:htgrid:int_lim)';
tgrid_rep = repmat(tgrid,1,I);
temp_muv = barw*params(5:I+4)./(params(5:I+4)+params(I+5:2*I+4)) + barw;
utvast_data = utildev_fun(temp_muv');
temp_tvast = utvast_data*tempM1(:,2:I+1);
compare = NaN*ones(T,ntgrid+1);
tt = 1; tempcheck = 1;
while tt <= T && tempcheck>=.5
    temp_tvast1 = repmat(temp_tvast(tt,:),ntgrid+1,1);
    compare(tt,:) = prod(((1+int_correc)*barw/(params(2)+2*params(1))< temp_tvast1 - tgrid_rep)...
        .*(temp_tvast1 - tgrid_rep<2*(1-int_correc)*barw/(params(2)+2*params(1))),2)';
    tempcheck = sum(compare(tt,:));
    tt = tt + 1;
end

if tempcheck<.5  || min(params(1:2*I+4))<.01^9 || params(2*I+7)<1
   out = -99^99;    lowu = NaN;
else
    % estimation of lowbaru
    lowbaru = min(utvast_data*[1; zeros(I,1)]);
    templl = NaN*ones(T,1);
    for t = 1:T
        temp5 = utvast_data(t,:);
        tempfun5 = @(x)cpdf_wv_integrand3(x,temp5(2:I+1)',temp5(1),I,params(2),...
            params(1),barw,params(5:I+4),params(I+5:2*I+4),params(2*I+5:2*I+6));
        intgrid = tgrid(logical(compare(t,2:(end-1))));
        temp6 = htgrid*(tempfun5(tgrid(1))+2*sum(tempfun5(intgrid))+tempfun5(tgrid(end)))/2;
        templl(t) = log(temp6) - (temp5(1,1)-params(3))^2/2/params(4) - log(params(4))/2 ...
            - log(1 - normcdf((lowbaru-params(3))/sqrt(params(4))));
    end
    
    out = sum(templl)/T + I*log(params(2)+2*params(1))+log(temp0c)-(I+1)*log(barw);
    lowu = lowbaru;
end
